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We discuss a representation of the Z3 Gauge-Higgs lattice field theory at finite density in terms of 
dual variables, i.e., loops of flux and surfaces. In the dual representation the complex action problem 
of the conventional formulation is resolved and Monte Carlo simulations at arbitrary chemical poten- 
tial become possible. A suitable algorithm based on plaquette occupation numbers and link-fluxes 
is introduced and we analyze the model at zero temperature and finite density both in the weak and 
strong coupling phases. We show that at zero temperature the model has different first order phase 
transitions as a function of the chemical potential both for the weak- and strong coupling phases. 
The exploratory study demonstrates that alternative degrees of freedom may successfully be used 
for Monte Carlo simulations in several systems with gauge and matter fields. 
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I. INTRODUCTORY REMARKS 

In the last three decades lattice QCD has seen an im- 
pressive development into a reliable tool for obtaining 
non-perturbative results in hadron physics. However, for 
one important application, the study of QCD at finite 
density, the lattice formulation has so far not lived up to 
our expectations. The reason is that for non-zero chem- 
ical potential the fermion determinant is complex and 
cannot be interpreted as a probability in a Monte Carlo 
simulation. Various approaches to circumvent this prob- 
lem were explored and the reviews at the annual lattice 
conferences provide a regular update [TJ |2] ■ 

An interesting approach to lattice systems with a com- 
plex action problem is to search for a mapping to alter- 
native variables where the partition sum is a sum over 
real and positive contributions such that in terms of the 
new variables the complex action problem is gone. Even 
for theories without a complex action problem an alter- 
native representation of the partition sum could allow for 
improved Monte Carlo simulations with new algorithmic 
ideas. A prominent simple example is the worm algo- 
rithm [3] for spin systems such as the Ising model which 
is based on a representation of the partition sum in terms 
of closed loops of Z2 flux. 

In recent years several examples for new algorithmic 
approaches that are based on transformations of lat- 
tice field theories to new variables were presented for a 
wide range of applications: Low dimensional systems [4], 
strongly coupled lattice field theories [5], systems with 
4-fermi interactions [5] , effective theories for QCD [3 [5] , 
scalar field theories 0HU] and U(l) lattice gauge theory 
112 ) . Similar in spirit, simulations directly based on 
the Trotter formula for lattice field theories in Hamilto- 
nian approach were explored [131 . Many interesting con- 
ceptual and algorithmic ideas emerged in these papers 
and systems that previously were not fully accessible to 
Monte Carlo simulations can now be explored. 



So far the studies with alternative variables (dual vari- 
ables) were mostly concerned with flux-like structures liv- 
ing on links, with the exception of the studies jTTl H2"] of 
pure U(l) gauge theory where the dual variables are sur- 
faces. In theories where gauge fields interact with matter 
fields, such as QCD, QED or Gauge-Higgs systems a dual 
representation will contain both: Surfaces for the gauge 
fields and fluxes for matter fields which serve as bound- 
aries for open surfaces. 

In this article we develop the idea of using dual repre- 
sentations for simulations of lattice field theories further 
and explore dual representations for the Z 3 Gauge-Higgs 
systems. This is a first step towards systems which couple 
gauge- and matter fields and where surfaces interacting 
with fluxes appear in the dual representation. As a mat- 
ter of fact all abelian Gauge-Higgs theories have a dual 
representation similar to the one we here discuss for the 
Z3 case and the techniques presented here can be adapted 
to other abelian cases. The choice to use the gauge group 
Z3 in this exploratory study is partly motivated by the 
possibility to couple a chemical potential and to explore 
finite density physics for Z3, which as the center group of 
SU(3) plays an interesting role in the phenomenology of 
QCD. At zero temperature we explore the various phase 
transitions as a function of the chemical potential. The 
main goal of this work is to develop further alternative 
representations of lattice field theories and their use in 
Monte Carlo simulations. 

In the next section we derive the dual representation 
of the partition sum in terms of fluxes and surfaces. In 
Section III we discuss observables and develop our strat- 
egy for the Monte Carlo simulation in terms of the dual 
variables. In Section IV we first compare the dual simu- 
lation of pure Z 3 lattice gauge theory and the full system 
at /1 = to results from a Monte Carlo calculation in 
the conventional approach. This is followed by the pre- 
sentation of the results at finite density (Section V). A 
summary in Section VI completes the paper. 
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II. DUAL REPRESENTATION OF THE Z 3 
GAUGE-HIGGS MODEL 

In this section we discuss the derivation of the dual rep- 
resentation for the Z3 Gauge Higgs model on the lattice. 
Both the Higgs- and the gauge variables are elements of 
the group Z 3 = {1, e i27r / 3 , e"" 1 '/ 3 }, and it is possible to 
couple a chemical potential /i. In the conventional repre- 
sentation the model then has a complex action problem 
at fi > 0. 

The Higgs field variables (p x live on the sites of a TV 3 x 
N t lattice and are parameterized as 4> x — e iSx2% ^ 3 with 
Sx € { — 1, 0, +1}. The gauge fields live on the links and 
are written as U Xi<7 = e ia ^,^/ 3 w ith a x>a E {-1,0, +1}. 
For both fields periodic boundary conditions are used 
in all 4 directions. The action S = Sg + Sh is split 
into a gauge action Sg and the action for the Higgs field 
Sh- For the gauge action we use Wilson's form (without 
constant term), 



S G 



+ u* 



(1) 



where the double sum runs over all plaquettes U Xit7T — 
U x , a U x+ -$, T U* +9 a U* T . The action for the Higgs field in 
a gauge field background configuration is 



>H 



(2) 

The hopping terms in the 4-direction (= temporal direc- 
tion) are coupled to the chemical potential /i. The factor 
r] controls the coupling between the Higgs- and the gauge 
fields. Essentially it plays the role of an inverse Higgs 
mass: If the Higgs field is infinitely heavy (77 = 0) it de- 
couples from the dynamics of the system. The partition 
sum is obtained by summing the Boltzmann factor over 
all possible configurations 



Z = 



£■ 

{s,a} 



-Sg—Sh 



(3) 



with 



z= n e 
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x s x — — 1 
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(4) 



For deriving the dual representation of the Higgs field 
action we use two identities which may be checked by 
explicit evaluation of both sides of the equations for the 
three cases s = — 1,0, 1. For the spatial nearest neighbor 
terms without chemical potential we use (s = —1,0, 1) 



exp(V^ s +ve^ s ) = C n 4 fe ' e ^ Sfe 
e 2r i + 2e~ r i 



fc=-i 



a 



e 2rt - e~ r ' 
e 2r i + 2e-" 



(5) 



For the temporal hops where the chemical potential en- 
ters we need (s = —1, 0, 1) 



exp 



2?7 cosh(ju) 



-f e -'T s 



£m 



fee' 3 



k=-i 



(6) 



+ 2e- r ' cosh ^ cos sinh(/i) -fey) 



These auxiliary formulas arc now used in an expansion 
of the Boltzmann factor. We find (the index j runs from 
1 to 3, V = iV 3 is the spatial volume and we use the 
representations 4> x = e ls * 27r /3 an( j jj x v _ e ia x ,„iir/3^ 

[Jexp (^e^^e 1 ^ 8 ^- 8 ^^ 



{fe} \x,j 



(7) 



We have introduced the shorthand notation J2{k} = 

Y\ x v Y^k x „=-i ^ or tne sum over au configurations of the 
expansion indices k. Inserting the expanded Boltzmann 
factor into the partition sum for the Higgs field in a given 
gauge field background (represented by the coefficients 

. X 

) we find 



Z H [a] =Y,z- Sh = Cf £ UbI^I )[l[M k 

W {k} \x,j 



n 
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where we have suitably reorganized the product over the 
link terms. The product of sums in the last term of this 
expression is the remaining sum over the configurations 
of the Higgs fields (parameterized by s x ). Each of the 
individual sums over the s x vanishes, unless [k x , v — 
k x -v,v\ is a multiple of 3. It is useful to define the triality 
function T(n) as 



Using the triality function T(n) we write the Higgs-field 



3 



partition sum as 



Z H [a 




(10) 



The last factor implements a constraint at each site x, 
i.e., the net-flux through x, given by [k XyL , — k x -v,u\ > 
has to vanish modulo 3. 

If the net flux through a site were to vanish exactly 
(not only modulo 3) the constraint would imply that only 
closed loops of flux are allowed. The fact that the net flux 
through each site vanishes only modulo 3 augments the 
closed loop condition with an additional rule allowing 3 
units of flux to emerge from a site or to vanish at a site. 

The next step is to use the flux representation (10) 



for the Higgs part in the expression for the full par- 
tition sum. The Boltzmann factor for the gauge part 
is also expanded, re-using the identity ([5| in the form 
(a = -1,0,1) 



exp 



C 



2 2 



e 15 + 2e- 
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p/3 - 
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Representing the Boltzmann factor for the gauge fields 
with this formula and inserting it in the full partition sum 
we find after some reordering of terms 

fn f n fn t ( e ^ - ^ j) 

xfn y^expu 2 ^^^ x 

\x.v a x u =—l \ 



_u<_cy 
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We have introduced plaquette occupation numbers p x>ar 
which may assume the values —1,0 and +1 and for unique 
labeling of the plaquettes we only consider p XyGT with 
a < t. As before, by J2{ P } we denote the sum over all 
configurations of the p x . aT . The last product of sums in 
(12) gives rise to triality constraints for all links which 



combine the plaquette occupation numbers of all plaque- 
ttes attached to that link and the fc-flux residing on that 
link. In a more compact notation the final result for the 
partition sum of the Z3 Gauge-Higgs model reads 



z = <7£5> 



[ P ,k}C F [k}W P [p}W F [k} . (13) 



The first sum runs over all configurations of the in- 
teger valued plaquette occupation variables p x ,ar € 
{ — 1, 0, +1} assigned to the plaquettes of the lattice, while 
the second sum is over all configurations of the flux vari- 
ables k X;V E {— 1, 0, +1} living on the links of the lattice. 
The flux variables k are subject to the constraint Cp[fc] 
given by 



C F {k) 



X \ V 



(14) 



which enforces the conservation of fc-flux modulo 3 at 
each site of the lattice (see ^ for the definition of T(n)). 
This flux conservation restricts the admissible configura- 
tions to closed oriented loops of fc-flux. A second con- 
straint, 



x,f v v<a 



Cp\p, k] = ^2[p x ,ua -Px-a.va] (15) 

^ \Px,au Px — a,ai^\ ~l~ k XyV j, 



connects the plaquette occupation numbers p with the fc- 
variables: At every link it enforces the combined flux of 
the plaquette occupation numbers attached to that link 
plus the fc-flux on that link to vanish modulo 3. Sim- 
ilar to the closed loop interpretation for the fc-flux, the 
link constraint forces the plaquette occupation numbers p 
into forming closed surfaces, or open surfaces bounded by 
loops of fc-flux, and again the surface rule is augmented 
by the modulo 3 exemption as for the fluxes. 

Both, the plaquette occupation numbers and the fluxes 
come with corresponding weight factors, given by 



w P \p}= n b 



|Px, CTT | 

p 



(16) 



x,a<r 



nv[fc]= ]Jb 




where the explicit expressions for B n and Bp are given 



in |5]) and (11). The partition sum (13) comes with an 
overall normalization factor C = (3 5 C^Ca) v where V = 



Nf 



and C v and Cp are given in 



51 and (11) 



III. OBSERVABLES AND MONTE CARLO 
UPDATE 

Having mapped the partition sums onto the dual vari- 
ables we now also need to identify the representation of 
the observables in terms of the dual degrees of freedom. 
In this exploratory study we concentrate on thermody- 
namical observables, i.e., observables that are obtained as 
derivatives of In Z with respect to the various couplings 
[20] . For the gauge sector we obtain the average plaque- 
tte (U) and the corresponding susceptibility \u a s first 
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and second derivatives with respect to the inverse gauge 
coupling ft, 



(U) 



1 9 1 7 

— In Z , XV 



1 



6NfN t d/3 



d 2 



InZ . (17) 



Both observables are normalized by the total number of 
plaquettes. The ^-dependence of the dual representation 



is encoded in the weight factors given by the product ( 16 ) 
over powers -BL P ' of the factors Bp defined in Jill) . There 
is an additional /^-dependence from the overall factor C 
which also contains Bp. This factor is necessary for the 
correct representation of the observables listed in (17 1. 

It is straightforward to evaluate the expressions (17) 
in the dual representation and one obtains: 



(U) = Bp + 




XU=B> 



6NfN t \B 



For the Higgs sector we study the particle number den- 
sity n and the corresponding susceptibility x-n which may 
be obtained as derivatives with respect to [i (again nor- 
malized with the 4- volume): 



d 



NfN t dfi 



InZ 



1 



d 2 



N*N t d/r 



InZ 



(19) 



The corresponding observables in the dual representation 
are 



1 



NfN t 
1 

NfNt 



(RiM + R Afo + R-iAf-i) , (20) 

(QiM + QoK + Q-iTVlr) 
R Af + i^A/li) 2 ) 



Ro^o 



- R-i 
R-i 



We introduced the abbreviations JV S , s = —1, 0, 1 for the 
total number of temporal link variables k x ^ that have a 
value s = —1,0,1. The R k and Q k with k — —1,0,1 
denote the ratios R k = M' k /M k and Q k = [M%M k - 
(M£) 2 ]/Mf , where the M k are the factors defined in M 
and the primes are used for their derivatives with respect 
to /x. 

Although the partition sums and the observables may 
seem somewhat involved in the dual representation, the 
dual Monte Carlo update turns out to be rather simple. 
We begin its discussion with introducing an update for 
pure gauge theory which in the next section we study as 



a first test case. The pure gauge update then serves as 
the starting point for developing the full algorithm where 
we will augment the pure gauge algorithm with another 
simple step to update the Higgs field as well. 

In the dual representation the pure gauge theory is a 
theory of only the plaquette variables and the partition 
sum can be written as (see (|13[)), 



(21) 



where the sum is over configurations of the plaquette 
occupation numbers p XMT G { — 1,0,1}. The constraint 



Cp[p] is simply the plaquette constraint (15) evaluated 



for the case when all fc-fluxes are set to k x 



0. The 



constraint is a product over all links and for each link 
the oriented flux of all plaquettes attached to that link 
has to be a multiple of 3 (see (15)). As discussed in the 
previous section the constraints give rise to the occupied 
plaquettes forming surfaces. Here, where we consider 
pure gauge theory, we have only closed surfaces since 
without matter fields there are no /c-fluxes that could 
serve as boundaries of open surfaces. 

For the update of the plaquette variables we start 
from a configuration of the plaquette occupation num- 
bers p x ,ar where all the link constraints are satisfied, us- 
ing the simplest choice, i.e., the trivial configuration with 
Px.ar = for all x and all tr, r. Starting from the triv- 
ial configuration we offer trial configurations which leave 
the constraints intact and accept them with the usual 
Metropolis probability. 

The trial configurations are generated by increasing 
or decreasing the plaquette occupation numbers on the 
faces of 3-cubes which we embed in four dimensions. Let 
1 < v\ < V2 < v % < 4 denote the three directions that 
define the 3-cube. Then there are four possible choices 
of the fi, i.e., four different 3-cubes can be embedded in 
four dimensions. Once the three directions Vi are fixed 
we select a site x for the lower left front corner of the 
cube and change the plaquette occupation numbers on 
the faces of the cube by ±1 according to one of the two 
possibilities depicted in Fig. [T] 

To take into account the fact that the plaquette oc- 
cupation numbers are restricted to —1, and +1, the 
addition of ±1 is understood only modulo 3. Addition 
modulo 3 is defined as the usual addition of the numbers 
+1,0 and —1, except for the two cases +1 + 1 = —1 and 
— 1 — 1 = +1. It is easy to check that the two changes 
illustrated in Fig. [T] leave the constraints at the sides of 
the cubes intact. The two possible changes are proposed 
with equal probability and are then accepted with the 
usual Metropolis probability. The corresponding accep- 
tance rate is a product of factors (Bp) ±x according to the 
plaquette occupation numbers at the faces of the cubes 
(see Eq. 16). A full cube sweep is a loop over all four 
possible embeddings of 3-cubes and over the N^ x N t 
possibilities to place a particular cube. 

On a finite lattice with periodic boundary conditions 
there are additional admissible configurations of the 



FIG. 1: Cube update: We illustrate the changes we propose 
for the plaquette occupation numbers at the faces of an em- 
bedded 3-cube with edges along the directions 1 < v\ < v 2 < 
vz < 4. The two choices in the lhs. and rhs. figure are pro- 
posed with equal probability. 



FIG. 2: Plaquette update: We illustrate the changes we pro- 
pose for the plaquette occupation numbers at the plaquette 
and the corresponding changes of the flux variables at the 
links of the plaquette. The directions 1 < v\ < v 2 < 4 de- 
termine the plane of the plaquette and the two choices in the 
lhs. and rhs. figure are proposed with equal probability. 



plaquette occupation numbers which are not properly 
treated by the cube updates alone. They consist of sur- 
faces that close around the periodic boundaries. To take 
them into account we offer an additional set of trial con- 
figurations, where we change all plaquette occupation 
numbers in a plane by ±1 (addition is again understood 
modulo 3) and accept this change with the correspond- 
ing Metropolis probability which is a product of pow- 
ers of Bp according to the plaquette occupation numbers 
on the plane of the trial configuration. For a complete 
" plane sweep" we offer this change for all possible planes 
(6 possible orientations of the planes with N% or N s N t 
slices for each orientation). For the update of pure gauge 
theory we mix plane and cube sweeps. 

For updating the full Gauge-Higgs systems also the 
matter fields need to be taken into account. They are 
represented by the link variables. The update of the k- 
fluxes is more interesting since they enter in two con- 
straints: At each site the total flux of the fc-variables has 
to vanish modulo 3. Furthermore, at each link the cor- 
responding fc-flux enters the constraint for the plaque- 
tte occupation numbers attached to that link. To deal 
with this situation we use proposal configurations that 
change a plaquette occupation number and the fc-ffuxes 
at the links of that plaquette. The two schoices for such a 
change are depicted in Fig. [2j where 1 < v\ < v% < 4 are 
the two directions that define the plane of the plaquette. 

It is easy to see that the two moves we propose do 
not change the constraints on the sites or on the links 
(addition is again only understood modulo 3 as for the 
cube and plane updates). The two changes in Fig. [2] are 
proposed with equal probability and are accepted with 
the Metropolis acceptance rate. This is a product of a 
ratio of the plaquette weights Bp and a product over 
ratios of the 4 link weights, i.e., powers of B n for the 
spatial links, and factors for the temporal links. A 
full plaquette sweep then consists of a loop over all 6 
embeddings of the planes of the plaquettes and all Nf x 
N t possibilities to place the plaquette. 

To summarize: For updating the full Gauge-Higgs 
model we use combined sweeps that mix full cube sweeps, 
full plane sweeps and full plaquette sweeps for the k- 
variables. 



It is interesting to note that when the Higgs field is 
coupled, i.e., for r\ > 0, the cube and plane sweeps could 
be omitted and it would be sufficient to work with the 
plaquette sweeps alone. The reason is that a cube up- 
date can be achieved by 6 plaquette updates, and a plane 
sweep by a combination of plaquette updates in the re- 
spective plane. However, we found that in particular 
omitting the cube updates may lead to a poorer perfor- 
mance of the algorithm in cases where the weights for the 
fc-fluxes are small. We remark, that for the Higgs model 
without gauge fields a very efficient worm algorithm [3 
can be constructed, and also for the full Z 3 -case a sur- 
face type of generalization of the worm idea is possible 
[T7] (see also [T^] for tests in this direction). 



IV. COMPARISON WITH CONVENTIONAL 
/i = RESULTS 

A. Pure gauge theory 

Having discussed the dual representation and the 
Monte Carlo algorithm in the dual picture we now come 
to the evaluation of the dual Monte Carlo approach. As 
a warm up exercise we begin with the case of pure gauge 
theory, where the Monte Carlo update is based on the 
cube and plane updates. The observables we consider are 
the plaquette expectation value (U) and the correspond- 
ing susceptibility xu- Their definitions and the corre- 



sponding dual expressions are given in (17) and (18 1. 

In Fig. [3] we show the results for (U) and xu in pure Z3 
lattice gauge theory (77 = 0), and compare the outcome 
of the dual simulation (circles) to the results from the 
conventional approach (crosses) using lattices with vol- 
umes of size 10 4 . The statistics is 10000 configurations 
separated by 10 cube sweeps and one plane sweep and an- 
other 10000 cube sweeps mixed with 1000 plane sweeps 
were used for equilibration. All errors we show in this 
work are statistical errors determined with the jackknife 
method. 

Fig. [3] shows that the results of the dual simulation 
perfectly agree with the outcome of the conventional ap- 
proach. Near /3 C ~ 0.7 the system apparently undergoes a 
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FIG. 3: Results for plaquette (U) (lhs. plot) and the plaquette susceptibility xu in pure Z3 gauge theory as a function of the 
inverse gauge coupling /3. We compare the conventional (crosses) and the dual approach (circles). 



quite prominent first order transition which separates the 
strong- {(i < (3 C ) and weak-coupling ((3 > (3 C ) phases, and 
also near the transition the results for the first and sec- 
ond derivatives of the free energy (i.e., (U) and \u) ob- 
tained with the conventional and dual approaches agree 
perfectly. We conclude that for the case of pure Z 3 lattice 
gauge theory the mapping to the dual representation, the 
identification of the observables and the simulation with 
the cube and plane algorithms work. 

It is interesting to inspect the acceptance rate for the 
plane sweeps. We find that the acceptance is zero for 
the strong coupling phase, i.e., for j3 < (3 C . For a small 
volume of size 4 4 we then see the onset of nonzero ac- 
ceptance at j3 c ~ 0.7 reaching a value of 0.28 at f3 = 1.0. 
Repeating the same analysis on volumes of size 10 4 we 
again find non-zero acceptance only above j3 c ~ 0.7, but 
the increase with /3 is much slower, and at j3 — 1.0 we still 
see an acceptance rate smaller than 0.01. We conclude 
that the non-trivially winding sheets of occupation are a 
finite size effect that very quickly dies out with increasing 
volume. 

At this point we remark, that the case of pure U(l) 
gauge theory in a dual representation which is similar to 
our Z 3 form [2T] has been studied numerically in [TTJ [T2] 
(see also the discussion in [2j ) . The algorithms in [11] are 
very similar to our updates (however, partly without the 
global plane updates). In [T^] defects (i.e., boundaries) 
are introduced for the surfaces and a generalization of 
the worm algorithm 3 to surfaces is explored. 



B. The full Z3 Gauge-Higgs model at 77 = 

Let us now come to the full Z 3 Gauge-Higgs model. 
Again we would like to test the dual approach and verify 



its validity by comparison with a conventional simulation, 
which is possible for fi = 0. 

For checking the correct implementation of the dual 
approach we compare the results for simulations at \i = 
using two different values of the coupling 77: 77 = 0.1 and 
77 = 0.5 and plot the observables as a function of j3. In 
Fig. [4] and Fig. [5] we show the results for (U) (lhs. plots) 
and xu (rhs.) as a function of (3. Fig. [I] is for 77 = 0.1 
and Fig. [5] for 77 = 0.5 and we compare the results from 
a conventional simulation (crosses) to the outcome from 
the dual approach (circles) and again use volumes of size 
10 and the same sequence and amount of the different 
update sweeps as for pure gauge theory. 

For 77 = 0.1 the first order transition persists and we see 
very little change in our observables (U) and xu when 
comparing the results to pure gauge theory. Again we 
observe that the results from the dual simulation and 
the conventional approach match very well. However, 
since at 77 = 0.1 and vanishing chemical potential /j the 
influence of the Higgs field seems to be small, we consider 
a second, larger value of 77. 

The results for (U) and xu at f] = 0.5 and [i = 0.0 
are shown in Fig. [5] We now observe quite a change 
in the behavior of the observables in comparison to the 
pure gauge and 77 = 0.1 cases. The phase transition 
has apparently disappeared and we only find a smooth 
crossover-type of behavior between the strong- and weak 
coupling phases. The maximum of the susceptibility xu 
has shifted to rather small values - the crossover takes 
place near (3 — 0.28. The important fact is, that also here 
at a larger value of 77, where obviously the Higgs field has 
a much stronger influence, the results from the conven- 
tional approach and the dual simulation agree very well, 
again confirming the correctness of the implementation 
of the dual approach. 
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FIG. 4: Results for (U) (lhs. plot) and \u in full Z3 Gauge-Higgs theory at j] = 0.1 and /1 = 0.0 as a function of the inverse 
gauge coupling /?. We compare the conventional (crosses) and the dual approach (circles). 
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FIG. 5: Same as in Fig(4] but now for rj = 0.5. 



V. THE Z 3 GAUGE-HIGGS MODEL AT FINITE 
DENSITY 

Let us now come to the more interesting case of finite 
density. Here conventional simulations fail and the full 
potential of the dual approach can be unveiled. Before 
we start with the presentation of Monte Carlo results 
we first discuss some characteristic features of the dual 
representation at finite density. 



A. Finite density dynamics in the dual 
representation 

The dual representation of the Z 3 Gauge-Higgs model 
uses two sets of degrees of freedom, the plaquette occu- 



pation numbers p and the fluxes k. For the analysis of 
the mechanisms that drive the systems at finite density 
it is useful to think a little bit about the dynamics of 
the dual variables, and this subsection is devoted to that 
task. 

The dual degrees of freedom assume values in 
{-1,0, +1}, i.e., p x ^ T e {-1,0, +1} and k X}U € 
{ — 1, 0, +1}. A trivial value of the plaquette occupation 
number, i.e., p Xl ar = comes with a Boltzmann factor 



of 1 (compare (16)), while non-trivial values p x .ar ± 1 



give rise to a factor of Bp < 1 (see ( 11 ) for the definition 
of Bp). Thus non-trivial values of plaquette occupation 
numbers p are suppressed by their Boltzmann factor. On 
the other hand configurations with many p x ^ ar ^ have 
a much higher entropy and (as always) the interplay of 
entropy and Boltzmann factor gives rise to the first order 
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transition of the pure gauge theory discussed in Subsec- 
tion IV-A. The corresponding observables (U) and xu 
are simple functions of the plaquette occupation numbers 
and their fluctuations. We stress at this point that both 
observables have a /3-dependent additive term (compare 
(18)). For the plaquette expectation value the additive 
term is given by Bp, and (U) is non- vanishing for /3 > 
even when all plaquette occupation numbers p are trivial, 
since Bp > for (3 > 0. 

Similar to the plaquette occupation numbers, the spa- 
tial flux variables k x j,j = 1, 2, 3, have a Boltzmann fac- 
and a Boltzmann factor B n < 1 



tor of 1 for fc Xi 
for k x ,j = ±1 (see (16)). As for the case of the pla- 
quette occupation numbers, we find for the fc-variables 
that trivial values of the spatial fluxes are preferred 
by the Boltzmann factor. The temporal flux variables 
k X! 4 are connected with the Boltzmann factors M s with 
s E {-1,0,1} defined in @. For p > we have 
M_|_i > M_i (see also the discussion below) and tempo- 
ral flux with k Xt i = +1 is favored over negative temporal 
flux, i.e., k x< 4 = — 1. 

To illustrate the physical picture in terms of the dual 
representation, in Fig. [6] we show a few low- lying excita- 
tions of the Z 3 Gauge-Higgs model in the dual representa- 
tion. Thick red lines oriented with arrows are used for the 
fc-flux and filled blue squares for non- vanishing plaquette 
occupation numbers, and the circles in the squares indi- 
cate the orientation of the plaquette according to the sign 
of the corresponding plaquette occupation number p x ^ aT . 
The simplest excitations (the lhs. diagram in Fig. [6]) are 
an occupied plaquette surrounded by flux. At each link 
the flux is compensated by the plaquette. Occupied pla- 
quettes with suitable relative orientation can be attached 
to each other (see the example in the center diagram). 
At the link where they are joined the flux is absent and 
the contributions from the plaquette occupation numbers 
compensate. Finally the excitation in the rhs. diagram 
makes use of the fact that flux and plaquette numbers 
need to vanish only modulo 3. Three units of flux emerge 
from a site, travel in time and then terminate at another 




'IQt 




FIG. 6: Examples of low- lying excitations in the dual repre- 
sentation of the Z3 Gauge-Higgs model. We use red thick lines 
with arrows for the k flux variables and blue squares (with a 
circle showing the orientation) for the plaquette occupation 
numbers p. The rules for admissible configurations dictate 
that at each site the total flux from the fc-variables has to be 
a multiple of 3. In addition for each link the combined flux of 
fc-variables and plaquette occupation numbers p also has to 
be a multiple of 3. 



site. The constraints are again saturated with plaquettes, 
such that at the central temporal link three plaquettes 
together obey the constraint. This excitation, which re- 
sembles a baryon, carries three factors of M+i and thus 
is enhanced by the chemical potential. 

Without chemical potential the Boltzmann factor for 
the excitations shown in Fig. [6] follows a simple rule: The 
more non- vanishing plaquette occupation numbers p and 
fluxes k a configuration has, the lower is the correspond- 
ing Boltzmann weight, and the dynamics discussed in 
Subsection IV-B. for the full Z 3 Gauge Higgs system at 
fj, = is again determined by the interplay of entropy 
and Boltzmann weight. 

The situation is more complex when the chemical po- 
tential p is turned on. Then the weights for the temporal 
flux variables k Xt 4 obey M + i > M_i and the probabil- 
ity for positive temporal fluxes is increased relative to 
the probability for negative fluxes. In the lhs. plots of 
Fig. [7] we show the ratios M+x/Mq and M_i/Mo as a 
function of p using ?y = 0.1. The Metropolis probabili- 
ties for accepting a step from fe X) 4 = to k Xt 4 — +1 and 
k X 4 — —1, respectively are given by max(l, M +1 /M ) 
and max(l, M_i/Mo). We see that for all /1 7^ we have 
M_|_i/Mo > M_i/Mq and temporal fluxes k Xl 4 = +1 are 
always favored over negative ones. For values of p up to 
p ~ 2.8 the discrepancy between M +1 /M and M_ 1 /M 
remains large and the non-zero chemical potential pumps 
positive temporal fc-flux into the system. However, this 
fc-flux has to be compensated by plaquettes which costs 
Boltzmann weight and dampens the increase of temporal 
flux with k Xt 4 = +1. This interplay between the pump- 
ing with p and the damping by the Boltzmann factor of 
the plaquette variables can give rise to a phase transition 
at some critical value of p, when the positive temporal 
fc-flux starts to dominate and drags along the plaque- 
tte occupation numbers. In the next subsection we will 
see that this is indeed the case for suitable values of the 
couplings. 

It is, however, important to note that for p > 2.8 
the two ratios M + i/Mq and M_i/Mq start to approach 
each other again and they both have limit 1 for large 
p. This implies that at large values of p the proba- 
bilities for temporal fluxes k X A = +1, k x ,4 = and 
fc X) 4 = —1 are equal, and the chemical potential does 
not favor an orientation for the k x ^. Thus for suffi- 
ciently large p the Boltzmann factor of the plaquette 
occupation numbers may again dominate the physics. 
This raises an interesting question that then has to be 
answered in this region: How can the particle number 
density n keep growing with p if the weights M s that 
drive the /i-dependence for small p become degenerate 
at large pi The answer lies in the flux representation of 
n given in (20): The numbers N s ,s € {+1,0,-1}, for 



temporal links with k x ^ — s are weighted with the fac- 
tors R s = M' s /M s . Also these ratios approach each other 
for large p and all three grow as e M (rhs. plot in Fig. J7|. 
Using N+i +A/q = N^N t we conclude from ( p0[ ) 

that for large p we have n oc R a oc r\e^ . In the plots for n 
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FIG. 7: Lhs.: The ratios M+i/Mo and M-i/Mq that determine the transition from temporal flux k Xl 4 — to k Xt 4 = +1 and 
k X: 4 — —1 as a function of /i at r\ = O.f . Rhs.: The coefficients R s = M' a /M s that determine the contributions of temporal flux 
k x ,4 = s with s £ { — 1, 0, +1} to the particle number density n in the dual representation ( 20 1. 



and Xn of Fig. 8 we display these limiting curves and find 
that the Monte Carlo data nicely approach the expected 
asymptotic behavior. We found similarly good asymp- 
totic behavior also for the weak coupling results shown 
in Fig. 11, but since there we use a smaller interval on 
the x-axis this is not entirely obvious from the plots. 

It is a remarkable feature of the dual representation 
that part of the behavior of observables (e.g., the asymp- 
totic behavior in the above example) is already encoded 
in the expansion factors of the partition sum and the rep- 
resentation of the observables in terms of dual variables. 

B. Finite density results at strong coupling 

Let us now come to the numerical results at finite den- 
sity. We numerically analyzed the finite density behavior 
for both values of rj that were considered in Section IV, 
i.e., r\ = 0.1 and rj = 0.5. For the rj = 0.5 case we 
did not find transitions as a function of ji and thus do 
not present the corresponding results in this exploratory 
study. Instead we focus on rj = 0.1 and study the sys- 
tem for two couplings on both sides of the transition lo- 
cated at /3 C ~ 0.7. We begin with the value (3 = 0.6 in 
the strong coupling region in this subsection, continuing 
with /3 = 0.8 in the weak coupling regime in the next 
subsection. 

For the finite density study we use lattices of sizes TV 3 x 
50, with N s ranging from 2 to 12. The reason for the large 
N t — 50 is that we want to study the system at zero 
temperature. Although with N t = 50 the temperature 
T = 1/50 = 0.02 in lattice units is not exactly zero, this 
value is much smaller than any other involved scale and 
constitutes a good approximation of T = 0. We consider 



the full set of observables discussed in Section III, i.e., in 
addition to the plaquette (U) and its susceptibility Xu 
we now also analyze the particle number density n and 
its susceptibility Xn- 

In Fig. |] we show the results for the four observables 
as a function of the chemical potential ji and compare 
runs on N% x 50 lattices with different spatial extents 
N s = 4, 8 and 12. The data for N s = 8 and N s = 
12 fall on top of each other and only the N s = 4 data 
show a slight discrepancy near /1 ~ 2.57 due to finite 
volume effects (the transition zone is shifted towards a 
slightly smaller ji). It is obvious that for fx — fi c ~ 2.57 
the system undergoes a first order transition: Both first 
derivatives, the plaquette expectation value (U) and the 
particle number density n show a clear discontinuity. The 
corresponding susceptibilities Xu and Xn diverge at /x c , 
which is, however, somewhat hard to see here since the 
transition is so sharp: As a matter of fact for very finely 
spaced values of /i near fj, c we find values for Xn that 
are three orders of magnitude larger than the scale used 
in the plot. A second indication that the transition is 
very narrow is the fact that only for N s — 4 we see small 
finite volume effects: The transition is shifted slightly to 
the left and appears somewhat rounded (at least for the 
first derivatives (U) and n). This trend towards visible 
finite size effects continues when using N s = 2 (data not 
shown) . 

Finally, a dataset at 12 3 x 100, which corresponds to 
an even lower temperature of T = 0.01 in lattice units, 
falls on top of the 12 3 x 50 data. We conclude that we 
reliably describe the situation at zero temperature, and 
that for 77 = 0.1 and /3 = 0.6 the system undergoes a very 
narrow first order transition at fi c ~ 2.57. 

From the fact that the plaquette undergoes such a dras- 
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FIG. 8: Results for ((7), \u, n and Xn in the strong coupling region of Z3 Gauge- Higgs theory (77 = 0.1 and /3 = 0.6) as a 
function of the chemical potential fj,. We compare the results for three different spatial volumes with N 3 — 4, 8 and 12. The 
dashed curves in the plots for n and \n are the asymptotic curves expected from the behavior of the coefficients M s for large 
fj, (see the discussion in Subsection V-A.). 



tic change at /i c we conclude that the gauge dynamics 
plays an important role in the transition. This is also re- 
flected in the finding that the critical value fi c depends on 
j3. For j3 = 0.6 we observed fi c ~ 2.57, while for f3 = 0.65 
it is at n c ~ 2.35, and fi c ~ 3.0 for j3 — 0.55. 

Finally we point out that for large fi the results for 
n and Xn approach the asymptotic curves oc r/ex.p(fi) 
(dashed lines in the plots) which we derived in Subsection 
V-A. from the behavior of the coefficients M s . 

To further understand the nature of the transition, in 
Fig.[9]we show 3-D illustrations of the plaquette and flux 
occupation numbers for a value of /it = 2.4 < fi c (top row 
of plots) and \x = 2.7 > /i c . We consider 3-dimensional 
sections of the 4-dimensional lattice with only spatial di- 
rections (lhs. pair of plots in each row) and 3-dimensional 
sections with the vertical direction being time (rhs. pair 



of plots). In each pair the lhs. section illustrates the pla- 
quette occupation numbers (plaquettes with p = +1 are 
blue, those with p = — 1 are red), while the rhs. section 
of each pair displays the corresponding fluxes (links with 
k = +1 are blue, those with k = — 1 are red). 

It is obvious that below fi c (top row of plots) the oc- 
cupation numbers for plaquettes and for fluxes are very 
small, while for fi > fi c (bottom) we see a large abun- 
dance of non-zero occupation numbers. The transition 
thus is between a dilute phase where all occupation num- 
bers are small and a condensed phase characterized by 
high occupation numbers for both, fluxes and plaquettes. 
Furthermore, a close inspection of the 3-D plot for the 
fluxes at /i > fj, c on the temporal section (lower right plot) 
shows that positive temporal flux (vertical blue lines in 
the plot) dominates, as expected from the finite density 
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FIG. 9: 3-D illustration of typical configurations of plaquette occupation numbers p and fluxes k in the strong coupling phase 
(P — 0.6,77 = 0.1) for /i = 2.4 < fj, c (top row of plots) and for /1 = 2.7 > /i c (bottom). We use 3-dimensional sections of 
the lattice embedded in 4 dimensions and show purely spatial sections (1st and 2nd plot in each row) and sections where the 
vertical direction is time (3rd and 4th plot). In each pair the lhs. section shows the non-trivial plaquette occupation numbers 
and we use blue for plaquettes with p = +1 and red for p = —1. Likewise, in the rhs. plot of each pair we show the non-trivial 
link variables with blue links for k = 1 and red for k = — 1. 



picture in terms of the dual variables discussed in the 
previous subsection. 

We remark at this point that the value of (U) ~ 0.35 
which we observe below /i c is mainly due to the constant 
term in the dual representation (18). Here we use j3 = 



0.6 for the inverse gauge coupling and the constant term 
has a value of 0.327. Thus also at very low plaquette 
occupation numbers p the plaquette has a value of (U) ~ 
0.35 at /? = 0.6. The jump of (U) ~ from 0.35 to almost 
1 at fi c is due to the condensation of non-trivial plaquette 
occupation numbers p as is obvious from Fig. [9j 

To further illuminate the picture of a condensation of 
the plaquette occupation numbers triggered by temporal 
/fi-fiux enhanced by we now look at plots for the oc- 
cupation numbers for plaquettes and fluxes as a function 



of the chemical potential fi. Fig. 10 shows the number of 
non-trivial spatial plaquettes P s , i.e., the total number of 
Px.ar 7^ 0, the number of temporal non-trivial plaquettes 
Pt {Px,a4: 7^ 0), the number of non-trivial spatial fluxes 
S {k x ,j 7^ 0), and the numbers of positive and negative 
temporal fluxes T + and T_ (k x ^ = +1 and k X 4 = —1). 
All these occupation numbers are normalized such that 
the maximally possible occupation number is 1. 

The occupation numbers for spatial (P s ) and temporal 
(Pt) plaquettes essentially vanish below fi c and then take 



a jump where both show an occupation of roughly 0.37, 
with the temporal plaquettes slightly enhanced. This be- 
havior clearly reflects the condensation of the plaquettes 
we had discussed above. Similarly the occupation num- 
ber S for the spatial flux nearly vanishes below fi c where 
it jumps to a finite value. The occupation numbers for 
positive and negative temporal flux, T + and T_ show an 
interesting behavior: Both essentially vanish below \x c . 
At the transition their degeneracy is lifted and both jump 
to different values. As the chemical potential is increased 
further the occupation numbers for positive and negative 
temporal flux both approach the value 1/3 (marked by 
a dashed horizontal line) which is the value one expects 
when the corresponding Boltzmann weights M s become 
degenerate (see the discussion in the previous subsection 
where we analyzed the behavior of the weight factors M s 
shown in Fig. [7]). 



C. Finite density results at weak coupling 

We continue with the discussion of the finite density 
behavior for the weak coupling region, i.e., for j3 = 0.8 
and again 7/ = 0.1. Here the situation is different since 
in the weak coupling phase already at fi = we have 
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FIG. 10: Occupation numbers for the dual variables at strong 
coupling (J3 = 0.6, 77 = 0.1) as a function of \x from our 12 4 x 50 
lattices. We show the number of non-trivial spatial plaque- 
ttes P s , i.e., the total number of p x , aT 7^ 0, the number of 
temporal non-trivial plaquettes Pt (p x ,a4 7^ 0), the number 
of non-trivial spatial fluxes S (k x j 7^ 0), and the numbers of 
positive and negative temporal fluxes T+ and T_ (k X) 4 = +1 
and k x ,4 = — 1, respectively). All occupation numbers are 
given as intensive quantities, normalized such that the max- 
imally possible occupation number is 1. The horizontal line 
marks the value 1/3 which the occupation numbers T + and 
T_ are expected to approach for large fi. 



(U) ~ 1, i.e., the plaquette occupation numbers are al- 
ready condensed. As in the previous subsection we con- 
sider the plaquette (U), its susceptibility xui the particle 
number density n and the corresponding susceptibility 
Xn- In Fig- [11] we show our results as a function of /i, 
comparing runs on x 50 lattices with different spatial 
extents N s = 4, 8 and 12. 

The susceptibility \n shows clearly that also in the 
weak coupling regime at /3 = 0.8 we find a phase tran- 
sition which is located at jit c ~ 1.35. Inspection of the 
particle number density shows that n has a discontinuity 
at (i c and we are again dealing with a first order tran- 
sition here. As for the transition in the strong coupling 
regime which we discussed in the last subsection, also 
here the finite volume effects are rather small and are 
visible only for the N s = 4 data. The transition is very 
narrow and as before we remark that the maxima in our 
data for Xn are much higher than the range used for the 
vertical axis. 

For the plaquette expectation value (U) the situation 
is different from the strong coupling phase of the last 
subsection. Since here we start with (U) ~ 1 below /x c , 
the change of (U) at fi c is rather unspectacular with the 
plaquette simply developing a mild slope. Correspond- 
ingly the susceptibility xu shows a slight drop above \x c 
as the plaquette numbers become completely saturated. 
We conclude that the transition in the weak coupling 
phase is predominantly driven by the matter fields, i.e., 



the flux variables k in the dual language. This finding 
is supported by the fact that changing the coupling f3 
from 0.8 to nearby values in the weak coupling regime 
(P = 0.75 and j3 — 0.85) has no noticeable effect on the 
value of \i c . This underlines the statement that the tran- 
sition is driven by the flux variables in the background 
of condensed plaquette occupation numbers p. 

Again we also compare the results for n and \n to 
the expected asymptotic behavior (dashed curves in the 
plot). We found that as in the strong coupling case the 
lattice data very nicely approach the expected asymp- 
totic behavior (although this is not very clearly visible 
for the range of /i- values chosen in Fig. 111. 



In order to analyze the nature of the transition in the 



weak coupling regime, in Fig. 12 we repeat the 3-D illus- 
tration of the plaquette occupation numbers and fluxes of 
the previous subsection. It is obvious that the plaquette 
occupation numbers are large on both sides of the tran- 
sition, and the fluxes undergo their transition in a con- 
densed medium of plaquette occupation numbers. The 
flux transition is again manifest in an abrupt increase 
of flux, although with a smaller amplitude than in the 
strong coupling phase. As before we see a dominance of 
positive temporal flux above fi c , i.e, vertical blue lines in 
the lower right plot. 

We conclude the discussion of the transition in terms 
of the dual variables by again analyzing the occupation 
numbers as a function of the chemical potential. In 
Fig. [13] we show the weak coupling results for non-trivial 
occupation numbers of spatial and temporal plaquettes 
(P s ,P t ), of spatial flux (S), and the occupation numbers 
for positive and negative temporal flux (T + , T_). We use 
the same definitions as in Fig. [I0]for the strong coupling 
phase. Obviously the plaquette occupation numbers are 
large for all values of and only a very mild change 
at [i c ~ 1.35 is visible. The flux variables show a small 
step at fi c which, however, is less pronounced than in the 
strong coupling case. Again we observe a strong splitting 
between the positive and negative temporal fluxes which 
reflects the influence of the chemical potential. At the 
largest values of \i we observe that the occupation num- 
bers for positive and negative temporal flux approach 
each other in agreement with the physical picture devel- 
oped on the discussion of Fig. [7] 



VI. SUMMARY AND DISCUSSION 

In this article we explored the possibility of a dual sim- 
ulation of gauge theories with matter fields. Although 
our study is for a simple model, the Z3 Gauge-Higgs 
model, it captures some of the features that are expected 
also for more interesting theories, in particular the ap- 
pearance of surfaces for the gauge fields and loops of flux 
for matter. In addition it was demonstrated that the 
complex action problem of the conventional representa- 
tion at non-zero chemical potential is solved in the dual 
approach. 
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FIG. 11: Results for (U), xu, n and \n m the weak coupling regime of Z3 Gauge-Higgs theory (77 = 0.1 and /3 = 0.8) as a 
function of the chemical potential fj,. We compare the results for three different spatial volumes with N s — 4, 8 and 12. The 
dashed curves in the bottom plots represent the asymptotic behavior of n and \n at large fi. 



A suitable Monte Carlo update was developed which 
properly treats the constraints of flux conservation at the 
sites and the surface constraints based at the links of the 
lattice (both constraints are modulo 3). In a detailed 
comparison in the pure gauge case and at vanishing chem- 
ical potential it was shown that the dual approach and 
the algorithm reproduce the results from a simulation 
in the conventional representation, thus establishing the 
validity of the dual approach. 

To test the approach at finite chemical potential /i 
where conventional techniques fail, we explored the be- 
havior of observables as a function of fi for two sets of 
couplings in the strong and weak coupling domains. In 
both cases we found first order transitions which were 
discussed not only based on usual observables, but also 
in terms of occupation numbers for the dual variables. 
We stress at this point that our two case studies in the 
strong and weak coupling regimes do of course not con- 



stitute a systematic analysis of the phase diagram for the 
Z3 Gauge-Higgs model - a task we do not aim at here. 

The techniques developed in this paper can easily be 
generalized to other Gauge-Higgs systems with abelian 
groups [18]. The dual degrees of freedom are again sur- 
faces and loops of flux, although the structure of the con- 
straints and the weight factors will differ for other abelian 
groups. Some of the aspects and properties found in the 
Z3 system might, however, turn out to be universal fea- 
tures of a dual approach: 

• The dual formulation represents the system using 
only gauge invariant degrees of freedom, i.e., suit- 
able occupation numbers for the plaquettes and the 
gauge invariant nearest neighbor terms of the mat- 
ter fields. 

• Part of the dynamics is encoded in expansion co- 
efficients of the partition sum and the observables 
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FIG. 12: 3-D illustration of typical configurations of plaquette occupation numbers p and fluxes k in the weak coupling phase 
(ft = 0.8, r\ — 0.1) for /i = 1.3 < (top row of plots) and for [t, = 1.5 > /i c (bottom). We use 3-dimensional sections through 
the lattice embedded in 4 dimensions and show purely spatial sections (1st and 2nd plot in each row) and sections where the 
vertical direction is time (3rd and 4th plot). In each pair the lhs. plot shows the non-trivial plaquette occupation numbers and 
we use blue for plaquettes with p — +1 and red for p = — 1. Likewise, in the rhs. plot of each pair we show the non-trivial link 
variables with blue links for k = 1 and red for k = — 1. 



(e.g., the asymptotic /i-behavior in the system stud- 
ied here). 

• Suitable Monte Carlo algorithms turn out to be 
rather simple, and at least for some cases a worm- 
type generalization to surfaces may be possible. 

• The dual representation is not only a tool to solve 
the complex action problem, but also allows for a 
conclusive discussion of the mechanisms at the vari- 
ous phase transitions in terms of the dual variables. 

It is obvious that the current results only present first 
steps towards the more important cases of non-abelian 
gauge fields or systems with fermions. Nevertheless we 
expect that some of the techniques developed in the dual 
approach to abelian Gauge-Higgs systems might prove 
useful also for these more interesting cases. 
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FIG. 13: Same as Fig. 11 now at j3 = 0.8, r/ = 0.1: We show 
the number of non-trivial spatial plaquettes P s , the number of 
temporal non-trivial plaquettes P t , the number of non-trivial 
spatial fluxes S and the numbers of positive and negative 
temporal fluxes T+ and Tl. The horizontal line marks the 
asymptotic value 1/3 for T+ and T_. 
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